#ifndef __PARAMETERS_H_
#include "tppmc/parameters.h"
#endif

#include <cmath>
#include <fstream>

#ifdef USING_LHAPDF
#include <LHAPDF/LHAPDF.h>
#else
#include "tppmc/Enumerations.h"
#endif

#include "tppmc/Particle.h"

using namespace std;

FitValues* VALUES::Values(int i) {
   if(mValues.find(i) not_eq mValues.end()) {
      return &(mValues[i]);
   } // if
   return NULL;
}

const FitValues* VALUES::Values(int i) const {
   std::map<int, FitValues>::const_iterator iter = mValues.find(i);
   if(iter not_eq mValues.end()) {
      return &(iter->second);
   } // if
   return NULL;
}

double FitValues::EvaluateTransversity(const tppmc::pythpart& parton) const {
   
   double x = parton.X();
   
   double f = N() *
   pow(x, Alpha()) *
   pow(1. - x, Beta()) *
   pow(Alpha() + Beta(), Alpha() + Beta()) /
   pow(Alpha(), Alpha()) /
   pow(Beta(), Beta());
   
   return f;
}

double FitValues::EvaluateSivers(const tppmc::pythpart& parton) const {
   double kperp = parton.Pt();
   // The form of the parameterisation is the same as for transversity
   return EvaluateTransversity(parton) *
   sqrt(2. * exp(1.)) *
   kperp / sqrt(M2()) * exp(-kperp * kperp / M2());
}


void VALUES::Populate() {
   mValues.insert(std::make_pair(LHAPDF::UP,
                                 FitValues(a_up, b_up, n_up, m2_up)));
   mValues.insert(std::make_pair(LHAPDF::DOWN,
                                 FitValues(a_down, b_down, n_down, m2_down)));
   mValues.insert(std::make_pair(LHAPDF::STRANGE,
                                 FitValues(a_strange, b_strange, n_strange, m2_strange)));
   mValues.insert(std::make_pair(LHAPDF::UBAR,
                                 FitValues(a_anti_up, b_anti_up, n_anti_up, m2_anti_up)));
   mValues.insert(std::make_pair(LHAPDF::DBAR,
                                 FitValues(a_anti_down, b_anti_down, n_anti_down, m2_anti_down)));
   mValues.insert(std::make_pair(LHAPDF::SBAR,
                                 FitValues(a_anti_strange, b_anti_strange, n_anti_strange, m2_anti_strange)));
//   mValues[LHAPDF::UP].Speak();
//   mValues[LHAPDF::DOWN].Speak();
//   mValues[LHAPDF::UBAR].Speak();
//   mValues[LHAPDF::DBAR].Speak();
//   mValues[LHAPDF::STRANGE].Speak();
//   mValues[LHAPDF::SBAR].Speak();
}

PARAMETERS::PARAMETERS() {
   // init all parameters
   kt2_average  = 0.25;
   ptq2_average = 0.2;
   
   Sivers.SetReadFile("tables/sivers_parameters.dat");
   Sivers.ReadParameters();
   
   Sivers.SetSetsReadFile("tables/sivers_sets.dat");
   Sivers.ReadSets();
   
   Transversity.SetReadFile("tables/transversity_parameters.dat");
   Transversity.ReadParameters();
   
   
   Collins.SetReadFile("tables/collins_parameters.dat");
   Collins.ReadParameters();
   
}


//========================================================= parameters
void SETS::ReadParameters(){
//   std::cout << "ReadParameters" << std::endl;
   ifstream in;
   in.open(FILE_READ); // fit 
   
   if(!in) {
      cerr << "PARAMETERS: Unable to open file with parameters! " << FILE_READ << std::endl;
   } else {
      //std::cout << "PARAMETERS: Open " << FILE_READ << std::endl;
   }
   
   in >>  parameters.a_up >>  errors.a_up;
   in >>  parameters.b_up >>  errors.b_up;
   in >>  parameters.n_up >>  errors.n_up;
   in >>  parameters.m2_up>>  errors.m2_up;
   
   in >>  parameters.a_down     >>  errors.a_down;
   in >>  parameters.b_down     >>  errors.b_down;
   in >>  parameters.n_down     >>  errors.n_down;
   in >>  parameters.m2_down    >>  errors.m2_down;
   
   in >>  parameters.a_anti_up >>  errors.a_anti_up;
   in >>  parameters.b_anti_up >>  errors.b_anti_up;
   in >>  parameters.n_anti_up >>  errors.n_anti_up;
   in >>  parameters.m2_anti_up>>  errors.m2_anti_up;
   
   in >>  parameters.a_anti_down     >>  errors.a_anti_down;
   in >>  parameters.b_anti_down     >>  errors.b_anti_down;
   in >>  parameters.n_anti_down     >>  errors.n_anti_down;
   in >>  parameters.m2_anti_down    >>  errors.m2_anti_down;
   
   in >>  parameters.a_strange >>  errors.a_strange;
   in >>  parameters.b_strange >>  errors.b_strange;
   in >>  parameters.n_strange >>  errors.n_strange;
   in >>  parameters.m2_strange>>  errors.m2_strange;
   
   in >>  parameters.a_anti_strange >>  errors.a_anti_strange;
   in >>  parameters.b_anti_strange >>  errors.b_anti_strange;
   in >>  parameters.n_anti_strange >>  errors.n_anti_strange;
   in >>  parameters.m2_anti_strange>>  errors.m2_anti_strange;
   
   
   in.close();
   
   main_parameters = parameters; // remember the main set of parameters
   main_parameters.Populate(); // Fill the mapped FitValues also
}


//========================================================= parameters
void SETS::UseSet(int k){
   parameters =  set[k];
}

//========================================================= parameters
void SETS::UseMainSet(void){
   parameters = main_parameters;
}


//========================================================= Constructor
void SETS::SetReadFile(const char* file){
   FILE_READ = file;
}

//========================================================= Constructor
void SETS::SetSetsReadFile(const char* file){
   SETS_FILE_READ = file;
}


//========================================================= Constructor
void SETS::SetWriteFile(const char* file){
   FILE_WRITE = file;
}

//========================================================= Constructor
void SETS::SetSetsWriteFile(const char* file){
   SETS_FILE_WRITE = file;
}

//========================================================= parameters
void SETS::ReadSets(){
   
   ifstream in;
   in.open(SETS_FILE_READ); // fit 
   
   if(!in) {
      cerr << "PARAMETERS: SETS: Unable to open file with parameters! " << SETS_FILE_READ << std::endl;
   } else {
      //std::cout << "PARAMETERS: SETS: Open " << SETS_FILE_READ << std::endl;
   }
   
   for(int i = 0; i < set_number; i++){
      in >>  set[i].a_up;
      in >>  set[i].b_up;
      in >>  set[i].n_up;
      in >>  set[i].m2_up;
      
      in >>  set[i].a_down;
      in >>  set[i].b_down;
      in >>  set[i].n_down;
      in >>  set[i].m2_down;
      
      in >>  set[i].a_anti_up;
      in >>  set[i].b_anti_up;
      in >>  set[i].n_anti_up;
      in >>  set[i].m2_anti_up;
      
      in >>  set[i].a_anti_down;
      in >>  set[i].b_anti_down;
      in >>  set[i].n_anti_down;
      in >>  set[i].m2_anti_down;
      
      
      in >>  set[i].a_strange;
      in >>  set[i].b_strange;
      in >>  set[i].n_strange;
      in >>  set[i].m2_strange;
      
      in >>  set[i].a_anti_strange;
      in >>  set[i].b_anti_strange;
      in >>  set[i].n_anti_strange;
      in >>  set[i].m2_anti_strange;
   }
   
   in.close();
}



//========================================================= parameters
void SETS::WriteParameters(){
   
   ofstream out;
   out.open(FILE_WRITE); // fit 
   
   if(!out) {
      cerr << "PARAMETERS: Unable to open file to write parameters! " << FILE_WRITE << std::endl;
   } else {
      std::cout << "PARAMETERS: Open " << FILE_WRITE << std::endl;
   }
   
   out <<  parameters.a_up << " " << errors.a_up << std::endl;
   out <<  parameters.b_up << " " <<  errors.b_up << std::endl;
   out <<  parameters.n_up << " " <<  errors.n_up << std::endl;
   out <<  parameters.m2_up<< " " <<  errors.m2_up<< std::endl;
   
   out <<  parameters.a_down     << " " <<  errors.a_down << std::endl;
   out <<  parameters.b_down     << " " <<  errors.b_down << std::endl;
   out <<  parameters.n_down     << " " <<  errors.n_down << std::endl;
   out <<  parameters.m2_down    << " " <<  errors.m2_down << std::endl;
   
   out <<  parameters.a_anti_up << " " <<  errors.a_anti_up << std::endl;
   out <<  parameters.b_anti_up << " " <<  errors.b_anti_up << std::endl;
   out <<  parameters.n_anti_up << " " <<  errors.n_anti_up << std::endl;
   out <<  parameters.m2_anti_up<< " " <<  errors.m2_anti_up << std::endl;
   
   out <<  parameters.a_anti_down     << " " <<  errors.a_anti_down << std::endl;
   out <<  parameters.b_anti_down     << " " <<  errors.b_anti_down << std::endl;
   out <<  parameters.n_anti_down     << " " <<  errors.n_anti_down << std::endl;
   out <<  parameters.m2_anti_down    << " " <<  errors.m2_anti_down << std::endl;
   
   
   out <<  parameters.a_strange << " " <<  errors.a_strange << std::endl;
   out <<  parameters.b_strange << " " <<  errors.b_strange << std::endl;
   out <<  parameters.n_strange << " " <<  errors.n_strange << std::endl;
   out <<  parameters.m2_strange<< " " <<  errors.m2_strange << std::endl;
   
   out <<  parameters.a_anti_strange << " " <<  errors.a_anti_strange << std::endl;
   out <<  parameters.b_anti_strange << " " <<  errors.b_anti_strange << std::endl;
   out <<  parameters.n_anti_strange << " " <<  errors.n_anti_strange << std::endl;
   out <<  parameters.m2_anti_strange<< " " <<  errors.m2_anti_strange << std::endl;
   
   
   out.close();
   
}

//========================================================= parameters
void SETS::WriteSets(){
   
   ofstream out;
   out.open(SETS_FILE_WRITE); // fit 
   
   if(!out) {
      cerr << "PARAMETERS: SETS: Unable to open file to write parameters! " << SETS_FILE_WRITE << std::endl;
   } else {
      std::cout << "PARAMETERS: SETS: Open " << SETS_FILE_WRITE << std::endl;
   }
   
   for(int i = 0; i < set_number; i++){
      out <<  set[i].a_up << " " ;
      out <<  set[i].b_up << " " ;
      out <<  set[i].n_up << " " ;
      out <<  set[i].m2_up << " " ;
      
      out <<  set[i].a_down << " " ;
      out <<  set[i].b_down << " " ;
      out <<  set[i].n_down << " " ;
      out <<  set[i].m2_down << " " ;
      
      out <<  set[i].a_anti_up << " " ;
      out <<  set[i].b_anti_up << " " ;
      out <<  set[i].n_anti_up << " " ;
      out <<  set[i].m2_anti_up << " " ;
      
      out <<  set[i].a_anti_down << " " ;
      out <<  set[i].b_anti_down << " " ;
      out <<  set[i].n_anti_down << " " ;
      out <<  set[i].m2_anti_down << " " ;
      
      
      out <<  set[i].a_strange << " " ;
      out <<  set[i].b_strange << " " ;
      out <<  set[i].n_strange << " " ;
      out <<  set[i].m2_strange << " " ;
      
      out <<  set[i].a_anti_strange << " " ;
      out <<  set[i].b_anti_strange << " " ;
      out <<  set[i].n_anti_strange << " " ;
      out <<  set[i].m2_anti_strange;
      out << std::endl;
   }
   
   out.close();
}


// Evaluate min and max of a function f(x) on the parameters sets
std::vector<double> SETS::MinMax( double (*f)(double ), double x ){
   vector<double> r;
   double min, max;
   
   //   for(int k = 0; k < set_number; k++){
   for(int k = 0; k < 2; k++){
      
      SETS::UseSet(k); // we do loop over all sets...
      
      double sechenie = (*f)( x );
      
//      std::cout << k << std::endl;
      
      if( k == 0 ){
         min = sechenie;
         max = sechenie;
      } else {
         if( (sechenie = (*f)( x ) ) > max ){
            max = sechenie;
         }
         if(  sechenie < min ){
            min = sechenie;
         }
      }
   }
   
   SETS::UseMainSet(); // main set is restored after the loop ...
   
   r.push_back(min);
   r.push_back(max);
   
   return r;
   
}

void PARAMETERS::Speak(std::ostream& os) const {
   os << std::endl; 
   os << "Transversity u-quark parameters:" << std::endl; 
   os << "N_u = " << Transversity.parameters.n_up << std::endl; 
   os << "alpha = " << Transversity.parameters.a_up << std::endl; 
   os << "beta = " << Transversity.parameters.b_up << std::endl; 
   os << "<k^2>_T = " << Transversity.parameters.m2_up << " GeV^2" << std::endl;
   os << std::endl; 
   os << "Transversity d-quark parameters:" << std::endl; 
   os << "N_d = " << Transversity.parameters.n_down << std::endl; 
   os << "alpha = " << Transversity.parameters.a_down << std::endl; 
   os << "beta = " << Transversity.parameters.b_down << std::endl; 
   os << "<k^2>_T = " << Transversity.parameters.m2_down << " GeV^2" << std::endl; 
   
   os << std::endl; 
   os << "Sivers u-quark parameters:" << std::endl; 
   os << "N_u = " << Sivers.parameters.n_up << std::endl; 
   os << "alpha = " << Sivers.parameters.a_up << std::endl; 
   os << "beta = " << Sivers.parameters.b_up << std::endl; 
   os << "M^2 = " << Sivers.parameters.m2_up << " GeV^2" << std::endl;
   os << std::endl; 
   os << "Sivers d-quark parameters:" << std::endl; 
   os << "N_d = " << Sivers.parameters.n_down << std::endl; 
   os << "alpha = " << Sivers.parameters.a_down << std::endl; 
   os << "beta = " << Sivers.parameters.b_down << std::endl; 
   os << "M^2 = " << Sivers.parameters.m2_down << " GeV^2" << std::endl; 
   os << std::endl; 
   os << "Sivers anti-u-quark parameters:" << std::endl; 
   os << "N_ubar = " << Sivers.parameters.n_anti_up << std::endl; 
   os << "alpha = " << Sivers.parameters.a_anti_up << std::endl; 
   os << "beta = " << Sivers.parameters.b_anti_up << std::endl; 
   os << "M^2 = " << Sivers.parameters.m2_anti_up << " GeV^2" << std::endl;
   os << std::endl; 
   os << "Sivers anti-d-quark parameters:" << std::endl; 
   os << "N_dbar = " << Sivers.parameters.n_anti_down << std::endl; 
   os << "alpha = " << Sivers.parameters.a_anti_down << std::endl; 
   os << "beta = " << Sivers.parameters.b_anti_down << std::endl; 
   os << "M^2 = " << Sivers.parameters.m2_anti_down << " GeV^2" << std::endl; 
   os << std::endl; 
   os << "Sivers s-quark parameters:" << std::endl; 
   os << "N_s = " << Sivers.parameters.n_strange << std::endl; 
   os << "alpha = " << Sivers.parameters.a_strange << std::endl; 
   os << "beta = " << Sivers.parameters.b_strange << std::endl; 
   os << "M^2 = " << Sivers.parameters.m2_strange << " GeV^2" << std::endl;
   os << std::endl; 
   os << "Sivers anti-s-quark parameters:" << std::endl; 
   os << "N_sbar = " << Sivers.parameters.n_anti_strange << std::endl; 
   os << "alpha = " << Sivers.parameters.a_anti_strange << std::endl; 
   os << "beta = " << Sivers.parameters.b_anti_strange << std::endl; 
   os << "M^2 = " << Sivers.parameters.m2_anti_strange << " GeV^2" << std::endl; 
   os << std::endl; 
   
   os << "Collins favored fragmentation parameters:" << std::endl; 
   os << "N_fav = " << Collins.parameters.n_up << std::endl; 
   os << "gamma = " << Collins.parameters.a_up << std::endl; 
   os << "delta = " << Collins.parameters.b_up << std::endl; 
   os << "M^2 = " << Collins.parameters.m2_up << " GeV^2" << std::endl;
   os << std::endl; 
   os << "Collins unfavored fragmentation parameters:" << std::endl; 
   os << "N_unfav = " << Collins.parameters.n_down << std::endl; 
   os << "gamma = " << Collins.parameters.a_down << std::endl; 
   os << "delta = " << Collins.parameters.b_down << std::endl; 
   os << "M^2 = " << Collins.parameters.m2_down << " GeV^2" << std::endl; 
   os << std::endl; 
}
